Install necessary packages: Set working directory

library(ggplot2)
library(readr)
library(sf)
Linking to GEOS 3.9.3, GDAL 3.5.2, PROJ 8.2.1; sf_use_s2() is TRUE
library(rnaturalearth)
The legacy packages maptools, rgdal, and rgeos, underpinning the sp package,
which was just loaded, will retire in October 2023.
Please refer to R-spatial evolution reports for details, especially
https://r-spatial.org/r/2023/05/15/evolution4.html.
It may be desirable to make the sf package available;
package maintainers should consider adding sf to Suggests:.
The sp package is now running under evolution status 2
     (status 2 uses the sf package in place of rgdal)
Support for Spatial objects (`sp`) will be deprecated in {rnaturalearth} and will be removed in a future release of the package. Please use `sf` objects with {rnaturalearth}. For example: `ne_download(returnclass = 'sf')`
library(rnaturalearthdata)

Attaching package: 'rnaturalearthdata'
The following object is masked from 'package:rnaturalearth':

    countries110
library(rgeos)
Loading required package: sp
rgeos version: 0.6-3, (SVN revision 696)
 GEOS runtime version: 3.9.3-CAPI-1.14.3 
 Please note that rgeos will be retired during October 2023,
plan transition to sf or terra functions using GEOS at your earliest convenience.
See https://r-spatial.org/r/2023/05/15/evolution4.html for details.
 GEOS using OverlayNG
 Linking to sp version: 2.0-0 
 Polygon checking: TRUE 
library(raster)
library(scales)

Attaching package: 'scales'
The following object is masked from 'package:readr':

    col_factor
library(tmap)

setwd('C:/Users/cy_su/PycharmProjects/DSCI_605_Data_Visualizations/Module 6/M6_Lab5/')

Lab 1: Load the Data:

plot_locations <- read_csv("HARV_PlotLocations.csv")
Rows: 21 Columns: 16
── Column specification ────────────────────────────────────────────────────────
Delimiter: ","
chr (11): geodeticDa, utmZone, plotID, stateProvi, county, domainName, domai...
dbl  (5): easting, northing, plotSize, elevation, plotdim_m

ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
raster_image <- raster("HARV_chmCrop.tif")

Lab 1: Convert raster to a data frame for ggplot2:

raster_df <- as.data.frame(rasterToPoints(raster_image))
names(raster_df) <- c("easting", "northing", "value")

Lab 1: Plot the Raster Image: Provide a legend title and better color gradient for the raster. Add an outline to the points and makes them more visible. Add a caption and refines labels. Apply a minimalistic theme for a cleaner look. Adjust text sizes and styles for better readability and aesthetics.

ggplot() +
  geom_raster(data = raster_df, aes(x = easting, y = northing, fill = value)) +
  scale_fill_gradientn(colors = terrain.colors(10), name = "Elevation") +
  geom_point(data = plot_locations, aes(x = easting, y = northing), color = "red", size = 4, shape = 21, fill = "white") +
  labs(title = "Plot Locations on Elevation Map", x = "Easting (m)", y = "Northing (m)", caption = "Source: HARV Dataset") +
  theme_minimal() +
  theme(
    plot.title = element_text(face = "bold", hjust = 0.5, size = 20),
    axis.title = element_text(face = "bold", size = 14),
    legend.title = element_text(face = "bold", size = 12),
    legend.text = element_text(size = 12),
    plot.caption = element_text(hjust = 0, size = 10)
  ) +
  coord_cartesian(expand = FALSE)

Lab 2: Load the Data:

watersheds <- st_read("Watersheds_HUC08_2009/WATERSHEDS_HUC08_2009_USDA_IN.shp")
Reading layer `WATERSHEDS_HUC08_2009_USDA_IN' from data source 
  `C:\Users\cy_su\PycharmProjects\DSCI_605_Data_Visualizations\Module 6\M6_Lab5\Watersheds_HUC08_2009\WATERSHEDS_HUC08_2009_USDA_IN.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 39 features and 12 fields
Geometry type: POLYGON
Dimension:     XY
Bounding box:  xmin: 374222.5 ymin: 4157935 xmax: 794277.8 ymax: 4691819
Projected CRS: NAD83 / UTM zone 16N
states <- st_read("tl_2019_us_state/tl_2019_us_state.shp")
Reading layer `tl_2019_us_state' from data source 
  `C:\Users\cy_su\PycharmProjects\DSCI_605_Data_Visualizations\Module 6\M6_Lab5\tl_2019_us_state\tl_2019_us_state.shp' 
  using driver `ESRI Shapefile'
Simple feature collection with 56 features and 14 fields
Geometry type: MULTIPOLYGON
Dimension:     XY
Bounding box:  xmin: -179.2311 ymin: -14.60181 xmax: 179.8597 ymax: 71.43979
Geodetic CRS:  NAD83

Set the theme:

tmap_mode("view")
tmap mode set to interactive viewing

Plot the Data: Set a color palette. Change the color of states to a simple grey to make the watersheds more visible For the watersheds adds a white border and sets up a pop-up showing the square miles of the area when an area is clicked on

palette <- RColorBrewer::brewer.pal(9, "YlGnBu")

tm_shape(states) +
  tm_polygons(col = "grey70", border.col = "grey40") +
  tm_shape(watersheds) +
  tm_polygons(col = "HUC_8", palette = palette, title = "Watershed (HUC_8)",
              border.col = "white", border.alpha = .5, id = "HUC_8",
              popup.vars = c("Area (sq miles)" = "SQ_MILES")) +
  tm_compass(type = "arrow", position = c("left","top"), size = 1.5) +
  tm_scale_bar(breaks = c(0, 50, 100), position = c("right","bottom"), size = 1) +
  tm_layout(main.title = "Watershed Map of Indiana (HUC_8)",
            main.title.position = "center",
            legend.position = c("left","bottom"),
            legend.text.size = 0.7,
            legend.title.size = 0.9,
            legend.bg.color = "white",
            legend.bg.alpha = 0.5,
            legend.frame = TRUE,
            frame = FALSE)
Warning: The argument size of tm_scale_bar is deprecated. It has been renamed to
text.size
Compass not supported in view mode.
legend.postion is used for plot mode. Use view.legend.position in tm_view to set the legend position in view mode.
Warning: Number of levels of the variable "HUC_8" is 39, which is
larger than max.categories (which is 30), so levels are combined. Set
tmap_options(max.categories = 39) in the layer function to show all levels.
Warning: In view mode, scale bar breaks are ignored.
LS0tDQp0aXRsZTogIk02X0xhYjVfU3BhdGlhbF9NYXBwaW5nX0NvZHlfWW9yayINCm91dHB1dDogaHRtbF9ub3RlYm9vaw0KLS0tDQoNCioqSW5zdGFsbCBuZWNlc3NhcnkgcGFja2FnZXM6KiogU2V0IHdvcmtpbmcgZGlyZWN0b3J5DQpgYGB7cn0NCmxpYnJhcnkoZ2dwbG90MikNCmxpYnJhcnkocmVhZHIpDQpsaWJyYXJ5KHNmKQ0KbGlicmFyeShybmF0dXJhbGVhcnRoKQ0KbGlicmFyeShybmF0dXJhbGVhcnRoZGF0YSkNCmxpYnJhcnkocmdlb3MpDQpsaWJyYXJ5KHJhc3RlcikNCmxpYnJhcnkoc2NhbGVzKQ0KbGlicmFyeSh0bWFwKQ0KDQpzZXR3ZCgnQzovVXNlcnMvY3lfc3UvUHljaGFybVByb2plY3RzL0RTQ0lfNjA1X0RhdGFfVmlzdWFsaXphdGlvbnMvTW9kdWxlIDYvTTZfTGFiNS8nKQ0KYGBgDQoNCioqTGFiIDE6IExvYWQgdGhlIERhdGE6KioNCmBgYHtyfQ0KcGxvdF9sb2NhdGlvbnMgPC0gcmVhZF9jc3YoIkhBUlZfUGxvdExvY2F0aW9ucy5jc3YiKQ0KcmFzdGVyX2ltYWdlIDwtIHJhc3RlcigiSEFSVl9jaG1Dcm9wLnRpZiIpDQpgYGANCg0KKipMYWIgMTogQ29udmVydCByYXN0ZXIgdG8gYSBkYXRhIGZyYW1lIGZvciBnZ3Bsb3QyOioqDQpgYGB7cn0NCnJhc3Rlcl9kZiA8LSBhcy5kYXRhLmZyYW1lKHJhc3RlclRvUG9pbnRzKHJhc3Rlcl9pbWFnZSkpDQpuYW1lcyhyYXN0ZXJfZGYpIDwtIGMoImVhc3RpbmciLCAibm9ydGhpbmciLCAidmFsdWUiKQ0KYGBgDQoNCioqTGFiIDE6IFBsb3QgdGhlIFJhc3RlciBJbWFnZToqKg0KUHJvdmlkZSBhIGxlZ2VuZCB0aXRsZSBhbmQgYmV0dGVyIGNvbG9yIGdyYWRpZW50IGZvciB0aGUgcmFzdGVyLg0KQWRkIGFuIG91dGxpbmUgdG8gdGhlIHBvaW50cyBhbmQgbWFrZXMgdGhlbSBtb3JlIHZpc2libGUuDQpBZGQgYSBjYXB0aW9uIGFuZCByZWZpbmVzIGxhYmVscy4NCkFwcGx5IGEgbWluaW1hbGlzdGljIHRoZW1lIGZvciBhIGNsZWFuZXIgbG9vay4NCkFkanVzdCB0ZXh0IHNpemVzIGFuZCBzdHlsZXMgZm9yIGJldHRlciByZWFkYWJpbGl0eSBhbmQgYWVzdGhldGljcy4NCmBgYHtyfQ0KZ2dwbG90KCkgKw0KICBnZW9tX3Jhc3RlcihkYXRhID0gcmFzdGVyX2RmLCBhZXMoeCA9IGVhc3RpbmcsIHkgPSBub3J0aGluZywgZmlsbCA9IHZhbHVlKSkgKw0KICBzY2FsZV9maWxsX2dyYWRpZW50bihjb2xvcnMgPSB0ZXJyYWluLmNvbG9ycygxMCksIG5hbWUgPSAiRWxldmF0aW9uIikgKw0KICBnZW9tX3BvaW50KGRhdGEgPSBwbG90X2xvY2F0aW9ucywgYWVzKHggPSBlYXN0aW5nLCB5ID0gbm9ydGhpbmcpLCBjb2xvciA9ICJyZWQiLCBzaXplID0gNCwgc2hhcGUgPSAyMSwgZmlsbCA9ICJ3aGl0ZSIpICsNCiAgbGFicyh0aXRsZSA9ICJQbG90IExvY2F0aW9ucyBvbiBFbGV2YXRpb24gTWFwIiwgeCA9ICJFYXN0aW5nIChtKSIsIHkgPSAiTm9ydGhpbmcgKG0pIiwgY2FwdGlvbiA9ICJTb3VyY2U6IEhBUlYgRGF0YXNldCIpICsNCiAgdGhlbWVfbWluaW1hbCgpICsNCiAgdGhlbWUoDQogICAgcGxvdC50aXRsZSA9IGVsZW1lbnRfdGV4dChmYWNlID0gImJvbGQiLCBoanVzdCA9IDAuNSwgc2l6ZSA9IDIwKSwNCiAgICBheGlzLnRpdGxlID0gZWxlbWVudF90ZXh0KGZhY2UgPSAiYm9sZCIsIHNpemUgPSAxNCksDQogICAgbGVnZW5kLnRpdGxlID0gZWxlbWVudF90ZXh0KGZhY2UgPSAiYm9sZCIsIHNpemUgPSAxMiksDQogICAgbGVnZW5kLnRleHQgPSBlbGVtZW50X3RleHQoc2l6ZSA9IDEyKSwNCiAgICBwbG90LmNhcHRpb24gPSBlbGVtZW50X3RleHQoaGp1c3QgPSAwLCBzaXplID0gMTApDQogICkgKw0KICBjb29yZF9jYXJ0ZXNpYW4oZXhwYW5kID0gRkFMU0UpDQpgYGANCg0KKipMYWIgMjogTG9hZCB0aGUgRGF0YToqKg0KYGBge3J9DQp3YXRlcnNoZWRzIDwtIHN0X3JlYWQoIldhdGVyc2hlZHNfSFVDMDhfMjAwOS9XQVRFUlNIRURTX0hVQzA4XzIwMDlfVVNEQV9JTi5zaHAiKQ0Kc3RhdGVzIDwtIHN0X3JlYWQoInRsXzIwMTlfdXNfc3RhdGUvdGxfMjAxOV91c19zdGF0ZS5zaHAiKQ0KYGBgDQoNCioqU2V0IHRoZSB0aGVtZToqKg0KYGBge3J9DQp0bWFwX21vZGUoInZpZXciKQ0KYGBgDQoNCioqUGxvdCB0aGUgRGF0YToqKg0KU2V0IGEgY29sb3IgcGFsZXR0ZS4NCkNoYW5nZSB0aGUgY29sb3Igb2Ygc3RhdGVzIHRvIGEgc2ltcGxlIGdyZXkgdG8gbWFrZSB0aGUgd2F0ZXJzaGVkcyBtb3JlIHZpc2libGUNCkZvciB0aGUgd2F0ZXJzaGVkcyBhZGRzIGEgd2hpdGUgYm9yZGVyIGFuZCBzZXRzIHVwIGEgcG9wLXVwIHNob3dpbmcgdGhlIHNxdWFyZSBtaWxlcyBvZiB0aGUgYXJlYSB3aGVuIGFuIGFyZWEgaXMgY2xpY2tlZCBvbg0KYGBge3J9DQpwYWxldHRlIDwtIFJDb2xvckJyZXdlcjo6YnJld2VyLnBhbCg5LCAiWWxHbkJ1IikNCg0KdG1fc2hhcGUoc3RhdGVzKSArDQogIHRtX3BvbHlnb25zKGNvbCA9ICJncmV5NzAiLCBib3JkZXIuY29sID0gImdyZXk0MCIpICsNCiAgdG1fc2hhcGUod2F0ZXJzaGVkcykgKw0KICB0bV9wb2x5Z29ucyhjb2wgPSAiSFVDXzgiLCBwYWxldHRlID0gcGFsZXR0ZSwgdGl0bGUgPSAiV2F0ZXJzaGVkIChIVUNfOCkiLA0KICAgICAgICAgICAgICBib3JkZXIuY29sID0gIndoaXRlIiwgYm9yZGVyLmFscGhhID0gLjUsIGlkID0gIkhVQ184IiwNCiAgICAgICAgICAgICAgcG9wdXAudmFycyA9IGMoIkFyZWEgKHNxIG1pbGVzKSIgPSAiU1FfTUlMRVMiKSkgKw0KICB0bV9jb21wYXNzKHR5cGUgPSAiYXJyb3ciLCBwb3NpdGlvbiA9IGMoImxlZnQiLCJ0b3AiKSwgc2l6ZSA9IDEuNSkgKw0KICB0bV9zY2FsZV9iYXIoYnJlYWtzID0gYygwLCA1MCwgMTAwKSwgcG9zaXRpb24gPSBjKCJyaWdodCIsImJvdHRvbSIpLCBzaXplID0gMSkgKw0KICB0bV9sYXlvdXQobWFpbi50aXRsZSA9ICJXYXRlcnNoZWQgTWFwIG9mIEluZGlhbmEgKEhVQ184KSIsDQogICAgICAgICAgICBtYWluLnRpdGxlLnBvc2l0aW9uID0gImNlbnRlciIsDQogICAgICAgICAgICBsZWdlbmQucG9zaXRpb24gPSBjKCJsZWZ0IiwiYm90dG9tIiksDQogICAgICAgICAgICBsZWdlbmQudGV4dC5zaXplID0gMC43LA0KICAgICAgICAgICAgbGVnZW5kLnRpdGxlLnNpemUgPSAwLjksDQogICAgICAgICAgICBsZWdlbmQuYmcuY29sb3IgPSAid2hpdGUiLA0KICAgICAgICAgICAgbGVnZW5kLmJnLmFscGhhID0gMC41LA0KICAgICAgICAgICAgbGVnZW5kLmZyYW1lID0gVFJVRSwNCiAgICAgICAgICAgIGZyYW1lID0gRkFMU0UpDQpgYGANCg==